\chapter{Rekonstruktionstheorie \index{Rekonstruktion} \index{Interpolation|see{Rekonstruktion}}\label{rekonstruktion}}
Die \textit{Rekonstruktion} - oder der hier synonym benutzte Begriff \textit{Interpolation} - bezeichnet Verfahren, die aus gegebenen diskreten Volumendaten die ursprünglichen (Mess-)Werte an einer beliebigen Position innerhalb des Volumens wiederherstellen. Nachfolgend wird speziell der eindimensionale Fall zur besseren Veranschaulichung besprochen. Die Abstraktion auf den dreidimensionalen diskreten Fall findet sich im Abschnitt \ref{anwendung}.

\section{Mathematische Grundlagen}
Zunächst benötigt es einige mathematische Begriffe, die die Basis der Rekonstruktionstheorie bilden. 
\subsection{Fouriertransformation\index{Fouriertransformation}\label{fourier}}
Eine \textit{Fouriertransformation} bildet eine integrierbare Funktion $g:\mathbb{R} \to \mathbb{R}$ auf ihre Fouriertransformierte $ \widehat{g} $, auch \textit{Spektrum} \index{Spektrum} von $g$ genannt, ab. Dabei gilt: \[\widehat{g}(x) := \frac{1}{\sqrt{2 \pi}} \int_\mathbb{R} g(t) e^{-itx} dt\] Die Inverse Fouriertransformation ist dann durch
\[g(t) = \frac{1}{\sqrt{2 \pi}} \int_\mathbb{R} \widehat{g}(x) e^{itx} dx\]gegeben. \cite[Abschnitt 2.1]{filter} \cite{wikiFourier}

\subsection{Faltung\index{Faltung}\label{faltung}}
Eine \textit{Faltung} zweier Funktionen $g,h: \mathbb{R} \to \mathbb{R}$ ist definiert durch \[(g \star h)(t) := \int_\mathbb{R} g(x)h(t-x) dx\] Dabei heißt $g$ die \textit{Signalfunktion} \index{Signalfunktion} oder kurz \textit{Signal} und $h$ heißt der \textit{Filterkern}, \textit{Rekonstruktionsfilter}, oder kurz \textit{Filter}. \index{Filterkern} \index{Filter} \cite[Seite 225]{real}

\subsection{Faltungssatz \index{Faltungssatz} \label{faltungssatz}}
Jede Fouriertransformation eines Produkts zweier Funktionen\footnote{Das Produkt zweier Funktionen $\varphi: \mathbb{R} \to A \subseteq \mathbb{R}$ und $\psi: \mathbb{R} \to B \subseteq \mathbb{R} $  ist definiert durch \\ $(\varphi\psi)(x) := \varphi(x) \cdot \psi(x) \quad \forall x \in \mathbb{R} $}  ist die Faltung der jeweiligen Fouriertransformation und umgekehrt:
\[ \widehat{gh} = \widehat{g} \star \widehat{h}, \quad \widehat{g \star h} = \widehat{g} \widehat{h} \]  \cite[Abschnitt 2.1]{filter}

\section{Rekonstruktion}
Mithilfe der eingeführten Definitionen und dem Faltungssatz ist es nun möglich, die Rekonstruktionstheorie aufzubauen: \\
Sei $g: \mathbb{R} \to \mathbb{R}$ ein ungesampeltes Signal, also eine kontinuierliche Funktion. Dann ist das Sampling äquivalent zu einer Multiplikation des Signals $g$ mit einem eindimensionalen Gitter, wie es Abbildung \ref{sampling} zeigt. Dazu definiert 
\begin{eqnarray*}
k:\mathbb{R} &\to& \lbrace 0, 1 \rbrace\\
t &\mapsto& \begin{cases} 1 \quad  t \in f\mathbb{Z} \\ 0 \quad \text{sonst} \end{cases}
\end{eqnarray*}
für eine Frequenz $f \in \mathbb{R}_{>0}$ ein solches Impulsgitter als Funktion. Diese ist eins bei allen Vielfachen von $f$ und null überall sonst. Das Produkt $gk$ identifiziert also das Sampling von $g$ mit der Samplingfrequenz $f$.
\vspace*{12pt}
\begin{figure}[bh]
\centering
\includegraphics[scale=1]{img/sampling.pdf}
\caption{Sampling}
\label{sampling}
\end{figure} \\
Die Fouriertransformation $\widehat{k}$ von $k$ ist wiederum eine Impulsfunktion mit der Frequenz $f$. Für die Fouriertransformation des gesampelten Signals $\widehat{gk}$ liefert der Faltungssatz die Gleichheit mit der Faltung von $\widehat{g}$ und $\widehat{k}$. Anschaulich wird daher $\widehat{g}$ an jeder Sampleposition - also mit der Frequenz $f$ - durch die Faltung mit $\widehat{k}$ wiederholt. (siehe Abbildung \ref{spektrum}) Die Kopie von $\widehat{g}$ im Zentrum um den Koordinatenursprung, also im Intervall $\left[-\frac{f}{2}, \frac{f}{2}\right]$, heißt \textit{Primärspektrum}\index{Primärspektrum}  von $\widehat{gk}$, alle anderen werden \textit{Aliasspektren} \index{Aliasspektrum}  von $\widehat{gk}$  genannt.
\begin{figure}[]
\centering
\includegraphics[scale=1]{img/spektrum.pdf}
\caption{Spektrum des Samplings}
\label{spektrum}
\end{figure} 
Ist $\widehat{g}$ dergestalt, dass sich Primär- und Aliasspektren von $\widehat{gk}$ nicht gegenseitig überlagern, kann aus $\widehat{gk}$ das ursprüngliche $\widehat{g}$ wiederhergestellt werden, indem die Aliasspektren so herausgefiltert werden, dass nur das Primärspektrum erhalten bleibt. Folglich ist eine Funktion $\widehat{h}$ gesucht, sodass die Multiplikation $\widehat{gk} \cdot \widehat{h}$ dieser Filterung entspricht.
\begin{figure}
\centering
\includegraphics[scale=1]{img/idealerFilterSpektrum.pdf}
\caption{Der ideale Filter im Spektralbereich}
\label{idealerFilterSpektrum}
\end{figure}
Ist $\widehat{h}$ - wie in Abbildung \ref{idealerFilterSpektrum} zu sehen - im Bereich des Primärspektrums eins und null sonst, erfüllt es trivialerweise diese Bedingung. Der Bereich wird dabei durch eine Frequenz $f_N < \frac{1}{2} f$, der sogenannten \textit{Nyquistfrequenz}, \index{Nyquistfrequenz} begrenzt. Damit entspricht $\widehat{h}$ einer Rechtecksfunktion mit dem Radius $f_N$ und es gilt \[\widehat{g} = \widehat{gk} \cdot \widehat{h}\]
Mithilfe der inversen Fouriertransformation und gleichzeitiger Anwendung des Faltungssatzes für die inverse Fouriertransformation auf der rechten Seite der Gleichung reduziert sich die Rekonstruktion eines gesampeltens Signals $gk$ folglich auf die Faltung mit dem Filterkern $h$: \[g = gk \star h\] \cite[Abschnitt 2.2]{filter}

\section{Idealer Filter\index{Idealer Filter}\label{ideal}}
Im Frequenzraum filtert die Rechtecksfunktion $\widehat{h}$, wie oben beschrieben, das Primärspektrum von den Aliaspektren von $\widehat{gk}$. Das herausgetrennte Primärspektrum entspricht der Fouriertransformation des Originalsignals. \\
Andererseits entspricht die Faltung $gk \star h$ genau der Rekonstruktion von $g$. Dementsprechend definiert die inverse Fouriertransformation der so gewählten Rechtecksfunktion $\widehat{h} = rect_{f_N}$ einen idealen Filter für $gk$.
\begin{eqnarray*}
h(t) &:=& \frac{1}{\sqrt{2 \pi}}\int_\mathbb{R} rect_{f_N}(x)e^{itx} dx
= \frac{1}{\sqrt{2 \pi}}\int_{-f_N}^{f_N} e^{itx} dx
= \frac{1}{\sqrt{2 \pi}} \left[\frac{e^{itx}}{it}\right]_{-f_N}^{f_N}\\[5pt]
&=& \frac{1}{i t \sqrt{2 \pi}} \left( e^{itf_N} - e^{-itf_N} \right)
= \frac{1}{i t \sqrt{2 \pi}} 2i\sin(tf_N)
= \frac{2}{\sqrt{2 \pi}} \frac{\sin(tf_N)}{t}\\[5pt]
&=& \sqrt{2 \pi} \frac{\sin(tf_N)}{t \pi}
\end{eqnarray*}
Für die Signalverarbeitung hat sich die Normierung \[ h(t) := sinc(\pi t) := \begin{cases} \frac{\sin(\pi t)}{\pi t} & \text{falls } t \neq 0 \\ 1 & \text{sonst} \end{cases} \] etabliert. \cite{wikiSinc} (siehe Abbildung \ref{sincFunktion} im Abschnitt \ref{abgSincFilter}) Da sich der Definitionsbereich von $sinc$ auf ganz $\mathbb{R}$ erstreckt, müsste in der Praxis für jede Rekonstruktion die Faltung auf dem gesamten Volumen durchgeführt werden, was selbst bei relativ kleinen Volumen nicht für eine Echtzeitanwendung praktikabel ist. Kapitel \ref{implementierung} behandelt daher Filter auf einem eingeschränkten Bereich.\\
Auch wenn der Begriff "`ideal"' etwas anderes suggeriert: Das Ergebnis ist, selbst bei einer Faltung über das gesamte Volumen, noch von der Beschaffenheit des gegebenen Volumens abhängig. Liegen im Spektrum des Volumens hohe Frequenzen vor, dann überlagern die Aliasspektren das Primärspektrum und es entsteht das sogenannte \textit{Aliasing}, \index{Aliasing} welches sich als Bildstörungen bemerkbar macht. \cite[Abschnitt 2.2, 2.3 und 3]{filter}

\section{Anwendung \label{anwendung}}
Für den praktischen Einsatz der Rekonstruktion eines Punktes $t$ eines Volumens reduziert sich das Integral der Faltung aus \ref{faltung} auf eine Summe, da das Signal $g$ typischerweise nur diskret vorliegt: \[(g \star h)(t) = \sum_{i=[t]-r}^{[t]+r} g(i)h(t-i) \] Aus Performancegründen wird die Faltung auf einen kleinen Radius $r$ um das nächstgelegene Voxel $[t]$ von $t$ eingeschränkt. In der Konsequenz werden die Filterfunktionen so gewählt, dass sie gerade außerhalb des Radius sehr klein oder null sind. Mit Ausnahme des klassischen Sinc Filtes, lässt sich dieser Sachverhalt leicht mithilfe der Funktionsplots der Filterkerne in Kapitel \ref{implementierung} einsehen. Durch diese Modifikation entstehen allerdings unausweichlich Bildfehler, sogenannte \textit{Artefakte}\index{Artefakt}.\\ 
Die grundlegende Vorgehensweise bei der Rekonstruktion eines Punktes $t$ ist es, alle umgebenen Voxel innerhalb des Radius zu bestimmen und mit dem Filter gewichtet aufzuaddieren. Abbildung \ref{rekonstruktionsSchema} illustriert das zugehörige Schema der Rekonstruktion mit dem Catmull-Rom Filter. Da der entsprechende Filterkern $h$ für $r > 2$ verschwindet, müssen in diesem Fall nur die Signalwerte $g_{-1}, g_0, g_1$ und $g_2$ für $g_i := g(i)$ mit jeweils $h(t - i)$ gewichtet, aufsummiert werden.\cite[Seite 225]{real} 
\begin{figure}[]
\centering
\includegraphics[scale=1]{img/rekonstruktion.pdf}
\caption{Rekonstruktion des Punktes $t$ mit einer Filterfunktion $h$}
\label{rekonstruktionsSchema}
\end{figure} \\
Alle behandelten Filter sind sogenannte \textit{separierbare} Filter, \index{Separierung} das heißt, die Erweiterung von ein- auf dreidimensionale Filter $h_r$ erklärt sich durch \[ h_r(x,y,z) = h_s(x) h_s(y) h_s(z)\] die Multiplikation der eindimensionalen Filter $h_s$ angewandt auf die einzelnen Komponenten des Punktes $(x,y,z)$. \cite[Abschnitt 2.3]{filter}